RungeKutta Function

private function RungeKutta(hini, Qin_begin, Qin_end, dt, geometry, rk, typOut, hTarget, eFlow, freeFlow, freeFlowElevation, QoutDOY) result(stage)

solve mass continuity of reservoir using runge kutta method

References:

Chow, V. T., Maidment, D. R. &Mays, L. W. (1988) Applied Hydrology. McGraw-Hill, New York, USA

Fenton, J.D. (1992) Reservoir routing, Hydrological Sciences Journal, 37(3), 233-246, DOI: 10.1080/02626669209492584

Arguments

Type IntentOptional Attributes Name
real(kind=float), intent(in) :: hini

initial stage value (m)
real(kind=float), intent(in) :: Qin_begin

input discharge at the beginning of time step (m3/s)
real(kind=float), intent(in) :: Qin_end

input discharge at the end of time step (m3/s)
integer, intent(in) :: dt

time step (s)
type(Table), intent(in) :: geometry

geometry of reservoir
integer, intent(in) :: rk

solution order 3 or 4
integer, intent(in) :: typOut

type ofoutflow: 1=free flow 2=target level
real(kind=float), intent(in) :: hTarget

follow a target (observed) stage (m)
real(kind=float), intent(in) :: eFlow

environmental flow (m3/s)
real(kind=float), intent(in) :: freeFlow

free flow (m3/s)
real(kind=float), intent(in) :: freeFlowElevation

free flow elevation (m asl)
character(len=3), intent(in) :: QoutDOY

doy to compute Qout

Return Value real(kind=float)


Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: Q1third
real(kind=float), public :: Q2third
real(kind=float), public :: Qout
real(kind=float), public :: area
real(kind=float), public :: dh
real(kind=float), public :: dh1
real(kind=float), public :: dh2
real(kind=float), public :: dh3
real(kind=float), public :: dh4
real(kind=float), public :: h1
real(kind=float), public :: h2
real(kind=float), public :: h3

Source Code

FUNCTION RungeKutta &
( hini, Qin_begin, Qin_end, dt, geometry, rk, typOut, hTarget, &
  eFlow, freeFlow, freeFlowElevation, QoutDOY ) &
RESULT (stage)


IMPLICIT NONE

! Function arguments
! Arguments with intent(in):
REAL (KIND = float), INTENT(IN)  :: hini !!initial stage value (m)
REAL (KIND = float), INTENT(IN)  :: Qin_begin !!input discharge at the beginning of time step (m3/s)
REAL (KIND = float), INTENT(IN)  :: Qin_end !!input discharge at the end of time step (m3/s)
INTEGER, INTENT(IN)              :: dt !!time step (s)
TYPE (Table), INTENT(IN)         :: geometry !!geometry of reservoir
INTEGER, INTENT(IN)              :: rk !! solution order 3 or 4
INTEGER, INTENT(IN)              :: typOut !!type ofoutflow: 1=free flow 2=target level
REAL (KIND = float), INTENT(IN)  :: hTarget !!follow a target (observed) stage (m)
REAL (KIND = float), INTENT(IN)  :: eFlow !!environmental flow (m3/s)
REAL (KIND = float), INTENT(IN)  :: freeFlow !! free flow  (m3/s)
REAL (KIND = float), INTENT(IN)  :: freeFlowElevation !! free flow elevation  (m asl)
CHARACTER (LEN = 3), INTENT(IN)  :: QoutDOY  !!doy to compute Qout 


! Arguments with intent(OUT):
REAL (KIND = float) :: stage

!local declaration:
REAL(KIND = float) :: h1, h2
REAL(KIND = float) :: h3  !used only if rk = 4
REAL(KIND = float) :: dh1, dh2, dh3
REAL(KIND = float) :: dh4	!used only if rk = 4
REAL(KIND = float) :: dh
REAL(KIND = float) :: Q1third, Q2third
REAL(KIND = float) :: Qout
REAL(KIND = float) :: area !(m2)


!------------end of declaration------------------------------------------------


SELECT CASE (rk)
   

   CASE(3)
        IF ( typOut == 2 .AND. hini < hTarget) THEN      
            
            Qout = MIN (eFlow, Qin_begin)
            
            !step 1
            !find area in table geometry
            CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', &
                                 keyOut='area', match='linear', &
                                 valueOut=area, bound = 'extendconstant' )
            
            !compute dh1
            dh1 = (Qin_begin - Qout) * dt / area
            
            !update h1
            h1 = hini + dh1 / 3.0
            
            !step 2
            Q1third = LinearInterp (0., REAL(dt), Qin_begin, Qin_end, 1.0/3.0*dt)
            
            !find area in table geometry
            CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', &
                                 keyOut='area', match='linear', &
                                 valueOut=area, bound = 'extendconstant' )
            
            !compute dh2
            dh2 = (Q1third - Qout) * dt / area

            !update h2
            h2 = hini + 2.0 * dh2 / 3.0

            !step 3
            Q2third = LinearInterp (0., real(dt), Qin_begin, Qin_end, 2.0/3.0*dt)
            
            !find area in table geometry
            CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', &
                                 keyOut='area', match='linear', &
                                 valueOut=area, bound = 'extendconstant' )
           
            !compute dh3
            dh3 = (Q2third - Qout) * dt / area
                 
        ELSE
          
            !step 1
            !find Qout 
            IF ( hini <= freeFlowElevation .AND. Qin_begin < freeFlow ) THEN
                Qout = Qin_begin
            ELSE
                CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', &
                                 keyOut = QoutDOY, match = 'linear', &
                                 valueOut = Qout, bound = 'extendconstant' )
                IF ( eFlow > 0.) THEN
                    Qout = MAX (eFlow, Qout)
                END IF
            END IF
            
           
            
            
            !find area in table geometry
            CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', &
                                 keyOut = 'area', match = 'linear', &
                                 valueOut = area, bound = 'extendconstant' )
            
            !compute dh1
            dh1 = (Qin_begin - Qout) * dt / area
            
            !update h1
            h1 = hini + dh1 / 3.0
  
            !step2
            Q1third = LinearInterp (0., REAL(dt), Qin_begin, Qin_end, 1.0/3.0*dt)
            
            !find Qout
            IF ( h1 <= freeFlowElevation .AND. Q1third < freeFlow ) THEN
                Qout = Q1third
            ELSE
                CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', &
                                 keyOut = QoutDOY, match = 'linear', &
                                 valueOut = Qout, bound = 'extendconstant' )
            END IF
            !find area in table geometry
            CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', &
                                 keyOut = 'area', match = 'linear', &
                                 valueOut = area, bound = 'extendconstant' )
            
            !compute dh2
            dh2 = (Q1third - Qout) * dt / area
            
            !update h2
            h2 = hini + 2.0 * dh2 / 3.0

            !step3
            Q2third = LinearInterp (0., real(dt), Qin_begin, Qin_end, 2.0/3.0*dt)
            
            !find Qout
            IF ( h2 <= freeFlowElevation .AND.  Q2third < freeFlow ) THEN
                Qout = Q2third
            ELSE
               CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', &
                                 keyOut = QoutDOY, match = 'linear', &
                                 valueOut = Qout, bound = 'extendconstant' )
            END IF
            !find area in table geometry
            CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', &
                                 keyOut = 'area', match = 'linear', &
                                 valueOut = area, bound = 'extendconstant' )
           
            !compute dh3
            dh3 = (Q2third - Qout) * dt / area
		END IF	  
  
        dh = dh1 / 4.0 + 3.0 * dh3 / 4.0

   CASE(4)
        !IF ( typOut == 2 .AND. hini < hTarget) THEN   
        !
        !    Qout = MIN (eFlow, Qin_begin)
        !    
        !    !step 1
        !    !find area in table geometry
        !    CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', &
        !                         keyOut='area', match='linear', &
        !                         valueOut=area, bound = 'extendlinear' )
        !    
        !    !compute dh1
        !    dh1 = (Qin_begin - Qout) *   area
        !    !write(*,*) dh1, Qout, eflow, Qin_begin, hini, htarget
        !
        !    !update h1
        !    h1 = hini + dh1 / 3.0
        !
        !    !step 2
        !    Q1third = LinearInterp (0., REAL(dt), Qin_begin, Qin_end, 1.0/3.0*dt)
        !    
        !    !find area in table geometry
        !    CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', &
        !                         keyOut='area', match='linear', &
        !                         valueOut=area, bound = 'extendlinear' )
        !    
        !    !compute dh2
        !    dh2 = (Q1third - Qout) * dt / area
        !    
        !    !update h2
        !    h2 = hini - 1.0 / 3.0 * dh1 + dh2
        !
        !    !step 3
        !    Q2third = LinearInterp (0., real(dt), Qin_begin, Qin_end, 2.0/3.0*dt)
        !    
        !    !find area in table geometry
        !    CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', &
        !                         keyOut='area', match='linear', &
        !                         valueOut=area, bound = 'extendlinear' )
        !    
        !    !compute dh3
        !    dh3 = (Q2third - Qout) * dt / area
        !
        !    !update h3
        !    h3 = hini + dh1 - dh2 + dh3
        !    
        !    !step4
        !    !find area in table geometry
        !    CALL TableGetValue ( valueIn = h3, tab = geometry, keyIn = 'h', &
        !                         keyOut='area', match='linear', &
        !                         valueOut=area, bound = 'extendlinear' )
        !
        !    !compute dh4
        !    dh4 = (Qin_end - Qout) * dt / area
        !    
        !ELSE
            
            !step 1
            !find Qout
            IF ( typOut == 2 .AND. hini < hTarget) THEN
                Qout = MIN (eFlow, Qin_begin)
            ELSE IF ( typOut == 2 .AND. hini >= hTarget) THEN
                CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', &
                                 keyOut = QoutDOY, match = 'linear', &
                                 valueOut = Qout, bound = 'extendconstant' )
                IF ( eFlow > 0. ) THEN
                    Qout = MAX (eFlow, Qout)
                END IF
                
            ELSE IF ( hini <= freeFlowElevation .AND. Qin_begin < freeFlow ) THEN
                Qout = Qin_begin
            ELSE
                CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', &
                                 keyOut = QoutDOY, match = 'linear', &
                                 valueOut = Qout, bound = 'extendconstant' )
            END IF
            
            !find area in table geometry
            CALL TableGetValue ( valueIn = hini, tab = geometry, keyIn = 'h', &
                                 keyOut = 'area', match = 'linear', &
                                 valueOut = area, bound = 'extendconstant' )
            
            !compute dh1
            dh1 = (Qin_begin - Qout) * dt / area
            
            !update h1
            h1 = hini + dh1 / 3.0
            
            !step2
            Q1third = LinearInterp (0., REAL(dt), Qin_begin, Qin_end, 1.0/3.0*dt)
            
            IF ( typOut == 1 ) THEN
            !update Qout
                IF ( h1 <= freeFlowElevation .AND. Q1third < freeFlow ) THEN
                    Qout = Q1third
                ELSE
                    CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', &
                                     keyOut = QoutDOY, match = 'linear', &
                                     valueOut = Qout, bound = 'extendconstant' )
                END IF
            END IF
            
            !find area in table geometry
            CALL TableGetValue ( valueIn = h1, tab = geometry, keyIn = 'h', &
                                 keyOut = 'area', match = 'linear', &
                                 valueOut = area, bound = 'extendconstant' )
       
            !compute dh2
            dh2 = (Q1third - Qout) * dt / area
            
            !update h2
            h2 = hini - 1.0 / 3.0 * dh1 + dh2
            
            !step3
            Q2third = LinearInterp (0., real(dt), Qin_begin, Qin_end, 2.0/3.0*dt)
            
            IF ( typOut == 1 ) THEN
            !update Qout
                IF ( h2 <= freeFlowElevation .AND.  Q2third < freeFlow ) THEN
                    Qout = Q2third
                ELSE
                   CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', &
                                     keyOut = QoutDOY, match = 'linear', &
                                     valueOut = Qout, bound = 'extendconstant' )
                END IF
            END IF
            !find area in table geometry
            CALL TableGetValue ( valueIn = h2, tab = geometry, keyIn = 'h', &
                                 keyOut = 'area', match = 'linear', &
                                 valueOut = area, bound = 'extendconstant' )
            
            !compute dh3
            dh3 = (Q2third - Qout) * dt / area
            
            !update h3
            h3 = hini + dh1 - dh2 + dh3
            
            !step 4
            
            IF ( typOut == 1 ) THEN
            !update Qout
                IF ( h3 <= freeFlowElevation .AND.  Qin_end < freeFlow ) THEN
                    Qout = Qin_end
                ELSE
                   CALL TableGetValue ( valueIn = h3, tab = geometry, keyIn = 'h', &
                                     keyOut = QoutDOY, match = 'linear', &
                                     valueOut = Qout, bound = 'extendconstant' )
                END IF
            END IF
            !find area in table geometry
            CALL TableGetValue ( valueIn = h3, tab = geometry, keyIn = 'h', &
                                 keyOut = 'area', match = 'linear', &
                                 valueOut = area, bound = 'extendconstant' )
           
            !compute dh4
            dh4 = (Qin_end - Qout) * dt / area
            
        !END IF
        
        dh = dh1 / 8.0 + 3.0 * dh2 / 8.0 + 3.0 * dh3 / 8.0 + dh4 / 8.0

   CASE DEFAULT
    
      CALL Catch ('error', 'RungeKutta', 'order Runge-Kutta not valid: ', &
	                      argument = ToString(rk))
	  STOP

END SELECT

!update reservoir stage
stage = hini + dh

RETURN
END FUNCTION RungeKutta